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Naval  Research  Laboratory  Ecological-Photochemical-Bio-Optical  Numerical 
Experiment  (Neptune)  Version  1:  A  Portable,  Flexible  Modeling  Environment  Designed 
to  Resolve  Time -Dependent  Feedbacks  between  Upper  Ocean  Ecology,  Photochemistry, 

and  Optics 


1.  Introduction 

For  over  98%  of  the  world’s  ocean  surface  area,  the  predominant  constituents 
contributing  to  surface  water  optical  properties  are  phytoplankton  and  non-living  organic 
matter  of  autochthonous  origin  [Morel  and  Prieur,  1977;  Morel,  1988;  Bricaud,  et  al., 
1998].  While  the  Case  I  (or  Bio-Optical)  assumption  presumes  that  these  constituents 
co-vary  in  space  and  time,  the  processes  that  constrain  the  abundance  of  phytoplankton 
and  non-living  organic  matter  are  fundamentally  different  such  that  optical  indices  for 
these  constituents  may  diverge  [ Carder ,  et  al.,  1989;  Siegel,  et  al.,  2002;  Siegel,  et  al., 
2005;  Lee  and  Hu,  2006],  Furthermore,  the  type  of  phytoplankton,  i.e.,  the  taxonomic 
composition  of  the  surface  phytoplankton  community,  may  further  impact  the  surface 
optical  property  variability  [e.g.,  Carder  and  Steward,  1985  \  Batch,  eta/.,  1996; 
Sathyendranath,  et  al.,  2004;  Westbury,  et  al.,  2005].  Specific  phytoplankton  groups  may 
be  adapted  to  specific  light,  nutrient,  and  hydrodynamic  regimes  and  their  dominance 
during  a  particular  season  or  within  a  particular  oceanographic  feature  may  further 
modulate  the  optical  properties  observed. 

Therein  lies  the  complexity  of  resolving  feedbacks  between  the  optical  properties  of  the 
surface  ocean  and  dynamic  changes  in  phytoplankton  abundance,  phytoplankton 
community  composition,  and  fluctuating  changes  in  the  optical  properties  and 
concentration  of  non-living  organic  matter:  there  is  a  continuing  time-dependent  synergy 
of  feedback  mechanisms  between  the  ecological  element  cycling,  which  directly  impacts 
the  outcome  of  phytoplankton  species  competition  as  well  as  overall  phytoplankton 
abundance,  and  processes  that  impact  the  optical  environment,  which  may  include  a 
combination  of  biochemical,  bio-optical,  photochemical,  and  hydrodynamic  mechanisms. 

The  NRL  Ecological-Photochemical-Bio-Optical  Numerical  Experiment,  or  Neptune  for 
brevity,  is  an  attempt  to  provide  a  flexible,  portable  modeling  construct  that  has  the 
capability  to  address  the  fundamental  optical  and  ecological  processes  that  detennine 
optical  property  variability  in  the  upper  ocean.  Neptune  is  fundamentally  based  upon  the 
established  lineage  of  ocean  biogeochemical  models  [e.g.,  Walsh,  1975;  Wroblewski, 
1977;  Walsh,  et  al.,  1978;  Jamart,  etal.,  1979;  Cullen,  1990;  Fasham,  et  al.,  1990;  Gregg 
and  Walsh,  1992;  Sarmiento,  et  al.,  1993;  Spitz,  et  al.,  2001]  wherein  transformations  of 
carbon,  nitrogen,  and  phosphorus  between  organic  and  inorganic  species  are  simulated  as 
a  consequence  of  lower  trophic-order  ecological  processes.  As  the  name  implies,  the 
modeling  construct  consists  of  three  main  components:  (1)  the  ecological  element  cycling 
computations  that  describe  phytoplankton  competition  for  and  transformations  of 
inorganic  nutrients  to  organic  matter  as  well  as  subsequent  organic  matter  cycling 
mediated  by  heterotrophic  organisms;  (2)  photochemical  transformations  of 
chromophoric  organic  matter  stimulated  by  incident  ultraviolet  irradiance;  and  (3)  a  bio- 
optical  module  that  computes  the  optical  properties  of  state  variables,  how  those  optical 
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properties  attenuate  spectral,  visible  light,  and  the  potential  photosynthetic  response  of 
photoautotrophic  organisms  to  the  absorption  of  a  budgeted  photon  flux  density. 

Flexibility  is  inherent  to  the  Neptune  modeling  construct;  no  a  priori  assumption  is  made 
that  a  single  version  of  the  model,  i.e.,  potential  material  flow  pathways  for  carbon, 
nitrogen,  and  phosphorus,  is  structurally  correct.  Rather,  multiple  modalities  of  material 
flow  may  be  examined  such  that  in  addition  to  parameter  space,  model  structural 
uncertainty  may  also  potentially  be  constrained  by  observational  fields.  Multiple 
functional  groups  of  photoautotrophs  and  heterotrophs  may  potentially  be  added  or 
reduced  to  a  single  composite  reservoir;  multiple  organic  matter  reservoirs  based  upon 
chemical  composition,  biodegradability,  or  optical  and  photochemical  properties  may  be 
added  or  restricted  as  well.  The  overall  philosophy  is  to  constrain  model  complexity 
down  to  an  essential  set  of  state  variables  and  parameters  requisite  to  provide  answers  to 
the  scientific  questions  being  addressed.  This  kind  of  flexibility  is  likely  mandatory  to 
begin  to  move  coupled  bio-optical/hydrodynamic  modeling  constructs  into  tactical  modes 
that  may  address  the  U.S.  Navy’s  desire  to  observe  and  predict  the  marine  environment  in 
support  of  U.S.  Navy  operations. 

1.1.  Core  Concepts 

Despite  this  flexibility,  Neptune  is  nonetheless  rooted  is  several  core  concepts  upon 
which  the  simulations  ultimately  rest.  First,  a  small  fraction  of  the  photons  absorbed  by 
non-living  organic  matter  and  phytoplankton  pigments  stimulate  photochemical  change: 
either  the  enzymatically  mediated  process  of  photosynthesis  or  various  abiotic 
photochemical  processes  [ Kieber ,  et  al.,  1989;  Mopper,  et  a/.,  1991].  In  either  case,  the 
second  law  of  photochemistry,  the  Stark-Einstein  Law,  states  that  a  single  photon  may 
only  potentially  stimulate  a  single  molecule  of  a  substance  to  react.  We  thus  reduce  both 
of  these  processes  down  to  a  simple  expression: 

(photons  absorbed)  ®  =  photochemical  yield  (1) 

where  <I>  is  the  apparent  quantum  yield  of  photosynthesis  or  the  apparent  quantum  yield 
of  abiotic  photochemical  reactions  involving  colored  dissolved  organic  matter  (CDOM). 
The  apparent  quantum  yield  is  simply  a  ratio  of  photons  absorbed  to  product  produced. 
Spectral  variability  is  approximated  by  summing  the  spectra  over  a  given  band  width,  Ak, 
of  discrete  increments: 


N 

photochemical  yield  =  ^  photons  absorbed(T)  ®(T)  AT  (2) 

i=l 

where  N  is  the  total  number  of  discrete  spectral  increments  over  a  given  spectral  range. 
Neptune  begins  at  the  full  ultraviolet  and  visible  spectral  range  of  290  -  700  nanometers, 
with  increments  beginning  at  10  nm  for  290  -  550  and  20  nm  for  550  -  700  nm. 

As  conservation  of  mass  and  energy  must  always  be  observed,  the  second  law  of 
photochemistry  allows  the  model  to  further  budget  how  photons  incident  upon  the 


Jolliff  and  Kindle 


Page  2  of  49 


ocean’s  surface  are  absorbed  over  discrete  depth  intervals  and  then  how  material 
reservoirs  within  those  intervals  are  subsequently  transformed.  In  the  case  of  CDOM, 
there  are  numerous  photochemical  transients  that  may  potentially  impact  the  biology  and 
ecology  of  the  upper  ocean  [Moran  and  Zepp ,  1997]  and  the  Neptune  construct  provides 
a  basic  scaffold  upon  which  those  processes  may  potentially  be  examined.  In  our  initial 
version  of  Neptune,  however,  we  are  concerned  primarily  with  the  photolytic  decay  of 
CDOM  because  the  reduction  of  CDOM  absorption  allows  additional  photons  to 
potentially  be  absorbed  by  phytoplankton  and  thus  stimulate  the  growth/accumulation  of 
organic  matter,  or  potentially  reduce  assimilation  rates  due  to  photoinhibition  of  the 
photosynthetic  process. 

Phytoplankton  assimilation  of  inorganic  carbon,  nitrogen,  and  phosphorus  does  not, 
however,  depend  solely  upon  the  absorption  of  photons,  but  may  be  constrained  by 
scarcity  of  other  requisite  resources.  Here,  as  in  many  other  ocean  biogeochemical 
models,  we  assume  that  Liebig’s  Law  of  the  Minimum  may  provide  a  conceptual 
construct  to  adequately  represent  phytoplankton  growth  based  upon  which  resource- 
based  growth  calculation  provides  the  lowest  organic  carbon  yield. 

Calculating  the  phytoplankton  carbon  assimilation  rate  based  upon  the  scarcity  of 
resources  requires  that  the  uptake  rate  of  nutrients  be  approximated  for  phytoplankton 
groups.  Uptake  of  nutrients  is  calculated  as  a  function  of  both  an  intrinsic  maximum 
assimilation  rate,  Vm,  and  Michaelis-Menten  uptake  kinetics  [Michaelis  and  Menten, 
1913]  where  affinity  for  nutrients  is  a  function  of  nutrient  concentration  and  a  half¬ 
saturation  constant,  Ks.  For  the  simple  case  of  the  change  in  product  concentration,  [P], 
as  a  function  of  substrate  concentration,  [S],  Michaelis-Menten  kinetics  are  described  by: 


d[P]  v  [S] 
dt  Ks  +[.9] 


(3) 


Another  core  assumption  is  that  phytoplankton  net  assimilation  of  carbon,  nitrogen,  and 
phosphorus  occurs  at  fixed  molar  ratios,  i.e.,  the  Redfield  [1963]  ratios  that  are  imprinted 
upon  the  deep  ocean.  Combining  Liebig’s  Law  of  the  Minimum,  Michaelis-Menten 
uptake  kinetics,  and  the  assumption  of  balanced  net  carbon  and  macronutrient 
assimilation  at  Redfield  proportions,  allows  one  to  replace  Vm  in  equation  (3)  with  a 
single  maximum  assimilation  rate,  pm,  and  calculate  an  array  of  potential  assimilation 
rates  wherein  each  array  element  is  potentially  reduced  by  the  scarcity  of  a  unique 
resource: 


ry  ] 

[Potential  Assimilation  Rates]  =  /um[ - ! - ],  i  =  1,  2,  3...  (4) 

Ks_i  '  [‘S'  ] 

The  minimum  of  the  array  is  then  the  realized  assimilation  rate,  in  accordance  with 
Liebig’s  Law  of  the  Minimum.  We  begin  by  describing  two  forms  of  nitrogen,  one 
reservoir  of  dissolved  inorganic  phosphorus,  and  a  spectrally  decomposed  photon  flux  as 
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potentially  limiting  agents,  but  this  could  be  potentially  expanded  to  include  dissolved 
forms  of  silica  and  micronutrients  such  as  iron. 

In  summary,  conservation  of  energy/mass,  the  second  law  of  photochemistry,  Leibig’s 
Law  of  the  Minimum,  Michaelis-Menten  uptake  kinetics,  and  Redficld  ratios  provide  the 
core  conceptual  framework.  The  link  between  upper  ocean  optics  and  ecology  is  then 
provided  by:  (1)  the  model  parameters  that  directly  relate  the  distribution  of  carbon 
within  various  organic  reservoirs  to  Inherent  Optical  Properties  (IOP’s);  (2)  the  equations 
that  describe  the  attenuation  and  scattering  of  photons  incident  upon  the  sea  surface  with 
depth;  and  (3)  the  equations  and  parameters  that  describe  how  absorption  of  photons 
stimulate  transformations  of  carbon,  i.e.,  flow  of  material  between  modeled  reservoirs. 

The  modeling  construct  presented  thus  unifies  ecological,  which  in  this  context  refers 
largely  to  the  competition  of  microorganisms  for  scarce  resources,  with  the 
photochemical/bio-optical  processes  into  a  single  simulation.  The  Neptune  modeling 
system  is  broken  down  into  several  components  that  are  coded  as  independent 
FORTRAN  90  subroutines.  The  purpose  of  this  report  is  to  document  these  subroutines, 
explain  the  science  behind  their  formulation,  and  provide  a  benchmark  for  later  reference 
as  the  model  is  modified  and  improved. 
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2.  Model  Details 


The  model  consists  of  a  main  program  unit  and  several  modules  containing  subroutines 
for  describing  both  boundary  conditions  and  various  ecological,  photochemical,  and 
optical  processes.  The  core  subroutine  functions  are  to:  (1)  define  the  spectral  irradiant 
boundary  conditions  for  visible  and  ultraviolet  bands  as  a  function  of  time-of-day,  day- 
of-year,  decimal  latitude,  and  decimal  longitude;  (2)  simulate  the  ecological  element 
transfonnations;  (3)  simulate  the  attenuation  of  spectral,  visible  irradiance  and  perform 
the  light-growth  phytoplankton  calculations;  (4)  simulate  the  attenuation  of  spectral, 
ultraviolet  irradiance  and  compute  the  abiotic  photochemical  yields  from  CDOM;  and  (5) 
provide  a  solution  to  the  simple  one-dimensional  case  where  state  variables  are  subjected 
to  vertical  diffusion  as  the  exclusive  hydrodynamic  process.  This  last  step  may  be 
discarded  when  Neptune  is  integrated  into  more  sophisticated  hydrodynamic  simulations. 
The  main  program  unit  is  responsible  for  declaring  and  initializing  the  state  variables  and 
some  key  parameters  that  are  sent  to  the  optical  and  ecological  subroutines.  The  main 
program  unit  then  updates  state  variables  after  the  ecological  and  photochemical  changes 
are  calculated  and  records  the  changes  using  netCDF  format  at  user  specified  intervals. 

The  number  of  state  variables  is  flexible  and  is  ultimately  based  on  how  many 
phytoplankton  functional  groups  and  classes  of  non-living  organic  matter  are  defined.  For 
the  initial  Neptune  start  up,  three  size-based  phytoplankton  functional  groups  and  three 
classes  of  non-living  organic  matter  are  defined.  A  size-based  phytoplankton  functional 
group  approach  was  selected  since  both  phytoplankton  optical  properties  and 
phytoplankton  affinity  for  nutrients  may  vary  with  gradations  in  overall  cell  size  [Raven, 
1998;  Bricaud,  et  al.,  2004]. 

With  the  exception  of  ammonium,  nitrate  (nitrite),  inorganic  phosphorus,  and  suspended 
sediments,  the  state  variables  are  tracked  in  units  of  carbon  concentration  and  the 
respective  nitrogen  and  phosphorus  contents  are  budgeted  by  using  fixed  C:N:P  ratios. 
With  the  notable  exception  of  heterotrophic  bacteria,  organic  reservoirs  become 
progressively  carbon  rich  subsequent  to  the  initial  organic  synthesis,  i.e.,  phytoplankton 
primary  production,  such  that  excess  nitrogen  and  phosphorus  budgeted  during  organic 
matter  transformations  are  returned  to  the  ammonium  and  inorganic  phosphorus 
reservoirs.  Since  heterotrophic  bacteria  are  nitrogen  and  phosphorus  enriched  compared 
to  their  organic  matter  substrates,  they  may  potentially  require  an  inorganic  nutrient 
supplement.  This  potential  requirement  is  dependent  upon  the  carbon  gross  growth 
efficiency  (GGE)  and  the  C:N:P  ratio  of  the  organic  substrate.  This  calculation  and  its 
structural  implications  for  the  simulated  ecosystem  are  discussed  in  Section  2.1. 

Figure  1  describes  the  base  Neptune  construct  for  state  variables  and  material  flow. 

Figure  2  shows  the  material  pathways  when  terrestrial/estuarine  boundary  fluxes  are 
included  with  additional  state  variables.  The  initial  material  flow  is  largely  dependent 
upon  the  grazing  fractions:  carbon  grazed  by  consumers  is  either  respired  or  diverted  to 
one  of  four  non-living  organic  matter  reservoirs.  Figure  3  shows  a  potential  structural 
modification  that  may  describe  colored  dissolved  organic  carbon  (CDOC)  production  as  a 
consequence  of  bacterial  processes. 
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Figure  1 .  Neptune- 1  material  flow  pathways:  circles  represent  processes,  rectangles  represent  state  variables, 
and  dashed  lines  are  remineralization  pathways. 
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Figure  2.  Neptune- 1  material  flow  pathways,  as  in  Figure  1,  with  added  terrestrial/fluvial/estuarine  boundary  fluxes. 
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and  dashed  lines  are  remineralization  pathways. 


In  the  first  construct  (Figures  1  and  2),  organic  matter  is  implicitly  produced  as 
semilabile.  The  second  (Figure  3)  shows  this  process  explicitly:  DOC  is  produced  at 
Redfield  [1963]  ratios  and  becomes  carbon  enriched  through  bacterial  modification.  This 
same  process  produces  some  CDOC.  We  show  these  two  model  structures  as  examples  of 
potential  modifications,  though  many  more  are  possible.  These  will  be  referred  to  as 
Neptune- 1  (Figure  1  and  2)  and  Neptune-2  (Figure  3). 

Below  we  describe  the  three  main  modeling  components:  ecological  element  cycling 
(2.1),  the  irradiant  boundary  condition  and  the  visible  optics  that  are  coupled  to  the 
phytoplankton  light-growth  calculations  (2.2),  and  the  photochemical  yield  subroutine 
(2.3).  The  remaining  section  describes  the  hydrodynamics  for  the  portable,  one¬ 
dimensional  case  that  is  coupled  to  the  MODAS  system  (2.4). 

2.1.  Ecological  Element  Cycling 

The  state  equation  for  each  state  variable  is  given  below,  followed  by  a  detailed 
description  of  the  ecological  computations.  The  main  programming  unit  for  describing 
element  cycling  and  phytoplankton  competition  for  nutrients  is  contained  in  the  module 
poseidon.f90.  The  specific  parameter  symbols  and  values  are  given  in  Tables  1,  2,  and  3. 
Parameter  values  are  based  on  literature  values  [Holm  and  Armstrong,  1981;  Ertel,  et  al., 
1986;  Fasham,  et  al.,  1990;  Benner,  et  al,  1992;  Chisholm,  1992;  Geider,  1992;  Kemp,  et 
al.,  1993;  Kirchman,  1994;  Vallino,  et  al.,  199 6;  Hopkins  on,  et  al.,  1998;  Walsh,  et  al., 
1999;  del  Giorgio  and  Cole,  2000;  Kirchman,  2000;  Mopper  and  Kieber,  2000;  Carlson, 
2002;  Hopkinson,  et  al.,  2002;  Wozniak,  et  al.,  2002;  Lucea,  et  al.,  2003;  Paytan,  et  al., 
2003],  unless  otherwise  noted.  Due  to  the  large  number  of  symbols  and  terms, 
supplemental  tables  that  describe  all  symbols  and  acronyms  may  be  found  in  the 
Appendix  to  aid  the  reader.  In  addition,  a  temperature  field  is  required  as  input  to  the 
model  since  growth,  grazing,  and  several  other  biological  functions  are  temperature 
dependent. 

All  phytoplankton  carbon  reservoirs  (P\,  i  =  1,  2,  or  3  in  Neptune  Version  1),  regardless 
of  how  many  phytoplankton  functional  groups  are  established,  conform  to  the  general 
equation: 


dt 


d  dP 

=  -V  •  {yp^)+—Kz—L+iinet  j 
oz  oz 


dWpbjPj 

dz 


(5) 


where  the  first  two  terms  on  the  right  hand  side  (RHS)  refer  to  the  hydrodynamic 
components.  The  first  term  (RHS)  is  the  advective  flux  divergence  (in  Cartesian 
coordinates),  which  is  shown  here  for  a  generic  tracer  X: 


V  •  (VX) 


r  du  dv 

- 1 - 1 - 

v  dx  dy  dz  y 


X 


(6) 
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Table  1.  Neptune- 1  Ecological  Model  General  Parameters  List 
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Table  2.  Neptune- 1  Phytoplankton  Functional  Group  Parameters 


Phytoplankton  Group  i=l,  Picophytoplankton  (0.2  -  2.0  um  cell  diameter)* 1 


gm3o  max  gross  carbon  assimilation  rate 

5.5  x10  V  @30°C 

n  |  MM  constant  phytoplankton  N03 

0.05  mmol  N  m"3 

n2  MM  constant  phytoplankton  NH4 

0.05  mmol  N  m"3 

n3  MM  constant  phytoplankton  DIP 

0.05  mmol  P  m'3 

wph3  maximum  phytoplankton  sinking  velocity 

O.Omd'1 

mimax  maximum  grazing  velocity 

2.3  x  10'5  s"1  @30°C 

Iv  Ivlev  grazing  paremeter 

0.5  m3  (mmol  C)"1 

a.1  fraction  of  grazed  P  carbon  respired 

0.6500 

a 2  fraction  of  grazed  P  carbon  to  DOC 

0.1890 

a3  fraction  of  grazed  P  carbon  to  SDET 

0.0792 

a4  fraction  of  grazed  P  carbon  to  LDET 

0.0008 

A 

a5  fraction  of  grazed  P  carbon  to  CDOC 

i-2>, 

»— i 

Phvtonlankton  Group  i=2.  Nanonhvtoplankton  12.0 

-  20.0  um  cell  diameter) 

gm3o  max  gross  carbon  assimilation  rate 

5.0  x10  V1  @30°C 

n  |  MM  constant  phytoplankton  N03 

0.17  mmol  N  m"3 

n2  MM  constant  phytoplankton  NH4 

0.17  mmol  N  m"3 

n3  MM  constant  phytoplankton  DIP 

0.17  mmol  P  m"3 

wphs  maximum  phytoplankton  sinking  velocity 

O.Omd1 

mimax  maximum  grazing  velocity 

2.3  x  10"5  s"1  @30°C 

Iv  Ivlev  grazing  paremeter 

0.3  m3  (mmol  C)"1 

a.i  fraction  of  grazed  P  carbon  respired 

0.5500 

a2  fraction  of  grazed  P  carbon  to  DOC 

0.1890 

a3  fraction  of  grazed  P  carbon  to  SDET 

0.0792 

a4  fraction  of  grazed  P  carbon  to  LDET 

0.0144 

A 

as  fraction  of  grazed  P  carbon  to  CDOC 

i-2>, 

t  —  r 

Phytoplankton  Group  i=3.  Microphvtoplankton  120.0  -  200.0  um  cell  diameter) 

gm3o  max  gross  carbon  assimilation  rate 

1.2  x  10'4  s"1  @30°C 

n  |  MM  constant  phytoplankton  N03 

1.3  mmol  C  m'3 

n2  MM  constant  phytoplankton  NH4 

0.9  mmol  C  m'3 

n3  MM  constant  phytoplankton  DIP 

0.5  mmol  P  m"3 

wpi,3  maximum  phytoplankton  sinking  velocity 

2.5  m  d'1 

nw  biomass  half-saturation  for  sinking  velocity 

0.75  mmol  C  m"3 

mimax  maximum  grazing  velocity 

2.3  x  10'5  s"1  @30°C 

Iv  Ivlev  grazing  paremeter 

0.3  m3  (mmol  C)"1 

ai  fraction  of  grazed  P  carbon  respired 

0.5000 

a2  fraction  of  grazed  P  carbon  to  DOC 

0.1540 

a3  fraction  of  grazed  P  carbon  to  SDET 

0.1400 

a4  fraction  of  grazed  P  carbon  to  LDET 

0.1400 

A 

a.5  fraction  of  grazed  P  carbon  to  CDOC 

i-2>, 

i=i 

*1  -  Parameter  selection  for  the  size-based  functional  groups  are  generally  within  the  range  of  published  values 
[e.g.,  Walsh  et  al.,  1999;  Fasham  et  al.,  1990],  but  calibrated  to  make  distinct  differences  in  accordance  with 
prevailing  theory  on  how  cell  size  and  surface -to-volume  ratios  impact  the  assimilation  of  and  affinity  for  nutrients 
as  well  as  intrinsic  rates  of  growth  [e.g.,  Chisholm,  1992;  Raven,  1998].  Grazing  fractions  are  based  on  the  DOM 
estimates  compiled  by  Carlson  [2002],  but  are  again  calibrated  so  that  large  phytoplankton  produce  more  particulate 
detritus,  and  proportionately  more  large  particulate  detritus,  than  do  smaller  phytoplankton.  The  Ivlev  parameter  is 
also  selected  based  on  size:  we  assume  that  consumers  of  smaller  phytoplankton  approach  maximum  consumption 
rates  at  lower  biomass  thresholds. 
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Phytoplankton  Group  i=l,  Picophytoplankton  (0.2 
a.1  fraction  of  grazed  P  carbon  respired 
02  fraction  of  grazed  P  carbon  to  labile  DOC 
03  fraction  of  grazed  P  carbon  to  SDET 
04  fraction  of  grazed  P  carbon  to  LDET 

Phytoplankton  Group  i=2,  Nanophytoplankton  (2,0 
Oi  fraction  of  grazed  P  carbon  respired 
02  fraction  of  grazed  P  carbon  to  labile  DOC 
03  fraction  of  grazed  P  carbon  to  SDET 
04  fraction  of  grazed  P  carbon  to  LDET 

Phytoplankton  Group  i=3,  Microphytoplankton 
Oi  fraction  of  grazed  P  carbon  respired 
02  fraction  of  grazed  P  carbon  to  labile  DOC 
03  fraction  of  grazed  P  carbon  to  SDET 
04  fraction  of  grazed  P  carbon  to  LDET 


0.1509  mol  N  (mol  C)_1 
0.0667  mol  N  (mol  C)'1 
0.0094  mol  P  (mol  C)_1 
-  mol  P  (mol  C)"1 

1.5  mmol  C  m'3 

49.5  mmol  C  m'3 

0.5 


General  Supplemental  Parameters 
03  i  molar  N  to  C  ratio  labile  DOM 
03  2  molar  N  to  C  ratio  semilabile  DOM 
C3  i  molar  P  to  C  ratio  labile  DOM 
C3  2  molar  P  to  C  ratio  semilabile  DOM 
n4  i  MM  constant  bacterial  labile  DOC 
n4  2  MM  constant  bacterial  semilabile  DOC 
Xcdoc  fraction  of  unassimilated  labile  DOC 
transferred  to  the  CDDC  reservoir 


(20,0  -  200.0  urn  cell  diameter) 
0.5000 
0.2200 
0.1400 
0.1400 


-  20.0  urn  cell  diameter) 
0.5500 
0.3564 
0.0792 
0.0144 


-  2.0  urn  cell  diameter) 
0.6500 
0.2700 
0.0792 
0.0008 
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The  horizontal  diffusivities  are  typically  ignored  due  to  the  large  inherent  numerical 
diffusion.  The  vertical  diffusivity,  Kz,  is  usually  a  prognostic  result  of  the  physical  model; 
however,  for  the  one-dimensional  case  Kz  may  be  approximated  from  temperature 
profiles  estimated  by  the  Modular  Ocean  Data  Assimilation  System  (MODAS).  This  is 
explained  in  section  2.4.  For  simplicity,  hydrodynamic  terms  are  dropped  for  the 
remaining  state  variable  equations. 

The  final  term  on  the  RHS  of  (5)  is  the  sinking  of  phytoplankton  cells;  the  sinking  rate  is 
most  often  set  to  0,  but  may  be  increased  to  parameterize  microphytoplankton  or  diatoms. 
The  third  term  on  the  RHS  of  (5)  is  the  net  phytoplankton  community  carbon  assimilation 
rate: 


=  (7) 

where  the  net  carbon  gain  for  phytoplankton,  junet,  is  the  maximal  realized  gross  carbon 
assimilation  rate,  gT,  minus  the  realized  basal  respiratory  rate,  s\,  and  realized  mortality 
tenn,  m\. 


The  maximal  realized  gross  carbon  assimilation  rate  for  phytoplankton  is  the  minimum  of 
the  calculated  light-limited  growth  rate,  gL,  the  nitrogen-limited  growth  rate,  gN,  and  the 
phosphorus-limited  growth  rate,  gp. 

8r  i  =min[gl  i’gN  i>gP  J  (8) 


Calculation  of  the  light-limited  (or  photo  inhibited)  growth  rate  is  discussed  in  Section  2.2 
and  coded  in  the  optics  module,  apollo.f90.  It  is  a  function  of  the  apparent  quantum  yield 
of  photosynthesis,  the  temperature-adjusted  maximal  gross  carbon  assimilation  rate,  and 
the  absorbed  photon  flux. 

gLJ  =  f(®c’grm.Tj’absorbed  _  photon  _  flux , )  (9) 


Absorption  and  maximal  temperature-adjusted  gross  carbon  assimilation  rates  are 
parameterized  differently  for  the  phytoplankton  functional  groups,  but  a  constant 
quantum  yield  is  used  for  all  groups.  This  may  be  changed  as  data  on  taxonomic  variation 
in  apparent  photosynthetic  quantum  yields  become  available. 


Nitrogen  limitation  is  a  function  of  nitrate  and  ammonium  uptake,  where  nitrate  uptake  is 
inhibited  in  the  presence  of  ammonium: 


g  V  i  g max  T  i 


[NO,] 


«i  ,.-hwo3] 


e-b[NHt]  + 


[NH4] 


n2  ,+[A//4] 


(10) 
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where  ti\  and  112  are  the  half-saturation  constants  for  nitrate  and  ammonium,  respectively, 
and  the  coefficient  b  is  the  nitrate  inhibition  factor  that  reduces  the  uptake  of  nitrate  in  the 
presence  of  ammonium. 


The  fraction  of  DIN  uptake  that  comes  from  nitrate  is  calculated  and  retained  for 
budgeting  purposes. 


X,  = 


TO 

vnLi^M93] 


'  [W03]  ' 

x  n\ ,-hwo,] 


e-b,lNH4]  + 


[ NH4 \ 


n2  ,+[AW4] 


4j; 


(ID 


Phosphorus  uptake  is  again  a  simple  Michaelis-Menten  function: 


SP  i  =  g  n 


[DIP] 


+[ DIP ] 


(12) 


The  maximum  assimilation  rate  at  30°  Celsius  is  adjusted  for  temperature  following  the 
expression: 


g- 


maxT  i 


=  min[gm30  , ,  ( gm30JeKtiJ-i0 ’ )]  (13) 


where  Kt  is  the  temperature  inhibition  coefficient  that  corresponds  to  a  Q10  value  of  1.88. 


A  basal  respiratory  cost  is  assigned  to  the  autotrophic  (and  hetero trophic)  biomass,  £\ 
(£2),  and  is  also  adjusted  for  temperature: 


(14) 


The  maximum  mortality  term  is  also  adjusted  for  temperature: 

^ImaxT =  minK max.i »  (  )l  (15) 

The  Ivlev  fonnulation  [Ivlev,  1955]  representing  a  single  consumer  group  utilizing  a 
single  food  source  determines  the  realized  grazing  rate: 

mu  =mimaxTj(l-e("(/V,)^))  (16) 

where  Ivi  is  the  Ivlev  parameter  that  determines  how  the  maximal  mortality  rate  is 
attenuated  with  decreasing  biomass.  The  concept  is  to  increase  grazing  stress  with 
increasing  phytoplankton  biomass  without  modeling  an  explicit  zooplankton  biomass. 
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Heterotrophic  bacterioplankton  conform  to  the  same  general  state  equation  as 
photoautotrophs,  but  without  any  potential  sinking  terms: 


dB_ 

dt 


=  Bne,bB 


(17) 


where 


Mnetb=(grb-£ 2  ~  mi)  (18) 

the  net  bacterial  growth  rate  is  the  maximal  realized  bacterial  assimilation  rate  minus  the 
basal  respiration  rate  and  the  mortality  rate.  Both  the  respiration  and  mortality  rates  are 
adjusted  for  temperature  as  in  phytoplankton  equations  (14)  and  (15).  The  model 
assumes  there  is  little  lag  between  bacterial  biomass  fluctuations  and  grazing  stress  from 
protozoan  consumers;  the  Ivlev  fonnulation  is  not  implemented: 

m2  =™2maxT(19) 

The  bacterial  maximum  assimilation  rate  is  again  the  function  of  the  least  available 
resource,  i.e.,  carbon,  nitrogen  or  phosphorus: 


Srb  nr  i  n  [  g  Ch ,  g  ^ ,  g  Pb  ]  (20) 

The  cellular  assimilation  of  DOC  is  a  product  of  the  temperature  adjusted  maximum 
gross  carbon  assimilation  rate,  the  carbon  gross  growth  efficiency  (GGE),  and  the  uptake 
kinetics  for  DOC  substrate: 


§Cb  =  gmaxT/,(GG£) 


"  [DOC]  N 
v«4  +[DOC]y 


(21) 


and 


gmaxT*  =min[gm3oA,(gm30AeS(^30))]  ^ 


To  calculate  the  assimilation  of  nitrogen  required  for  balanced  bacterial  growth,  a  ratio  is 
constructed  of  the  maximum  amount  of  inorganic  nutrients  that  could  potentially  be 
assimilated  by  the  heterotrophic  biomass  and  the  amount  of  NH4  supplement  required  if 
dissolved  organic  nitrogen  (DON)  is  insufficient  to  meet  balanced  growth  needs: 


§Nb  §Cb 


f  Maximum  -  Bacterial  -  NHA  -  Assimilation{MBNA )  ' 
Bacterial  -  NH4  -  Required  ( BNR ) 


(23) 
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where 


MBNA  = 


[NH4] 

U5+[a^4]J 


gm axX^^(24) 


BNR  = 


'  [DOC]  ' 
KnA+[DOC]y 


gmaxTh[GGEd2-0}]B(25) 


If  dissolved  organic  nitrogen  uptake  is  sufficient  to  meet  balanced  growth  needs  for 
bacteria  (C:N  =  5:1,  molar),  the  calculation  is  discarded  and  gm  =  gcb ■  At  a  GGE  of  30% 
and  a  semilabile  organic  matter  C:N  ratio  of  15,  the  calculation  is  not  performed  as 
organic  nitrogen  uptake  is  presumably  sufficient.  The  model  retains  the  flexibility  for 
inorganic  nitrogen  supplement  by  bacteria  should  these  parameters  change.  A  discussion 
of  this  type  of  calculation  may  be  found  in  Fasham  et  al.  [1990]. 


A  similar  calculation  is  perfonned  for  the  potential  inorganic  phosphorus  supplement. 
This  calculation  is  seldom  discarded  as  DOM  tends  towards  low  phosphorus  content  (C:P 
250:1)  and  bacterial  biomass  tends  towards  a  high  phosphorus  content  (C:P  -53:1). 


gpb  g  Cb 


f  Maximum  -  Bacterial  -  DIP  -  Assimilationi  M B PA )  ' 
Bacterial  -  DIP  -  Required(BPR) 


(26) 


MB  PA  = 


[DIP] 

n5  +  [DIP] 


V"5 


A  max  Tb  B%2 


(27) 


BPR  = 


[DOC] 
v«4  +  [DOC]y 


gM(GGE)B42 


(28) 


The  net  carbon  assimilation  is  then  calculated  according  to  equation  (20).  The  net  result 
of  these  calculations  is  that  the  heterotrophic  biomass  competes  with  the  autotrophic 
biomass  for  phosphate  when  DOC  is  abundant,  a  pattern  consistent  with  the  observations 
that  DOC-rich  aquatic  systems  tend  towards  phosphorus-limitation  [Kirchman,  1994]. 

For  Neptune- 1,  a  single  DOC  source  is  the  carbon  substrate  for  bacterial  growth. 
Inorganic  nutrient  supplement  may  be  required  depending  upon  how  the  C:N  and  C:P 
ratios  for  DOC  are  parameterized.  Neptune-2  assumes  a  highly  labile,  near-Redfield 
source  of  carbon  substrate  that  does  not  require  inorganic  mineral  nutrient  supplement 
(DOCi)  and  a  secondary,  semilabile  source  that  may  require  a  supplement  for  balanced 
bacterial  biomass  growth  (DOC2).  The  Neptune -2  bacterial  growth  calculations  are  then 
modified  to  account  for  two  potential  carbon  sources: 
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So,  =  min  [gmaxT/l  (GGE)DOCx,  gmaxTA  (GGE)] , 


DOCx  = 


(  [DOC,]  " 

_|_ 

r  [DOC2] 

[n^+iDOC,]  j 

i 

{ n42+[DOC2\ j 

(29) 


and  the  DOC  uptake  fractions  must  also  be  retained  for  proper  budgeting: 


Y  — 
^  DOC 


l DOC ,] 
n4_  \  +[DOCl\J 


l DOC ,] 
n4j  +[DOCJ  j 


[DOC] 
n42  +  [DOC2  ]  j 


T  (30) 


The  BNR  and  BPR  computations  are  then  also  adjusted  to  account  for  two  different 
sources  of  nitrogen  and  phosphorus: 


LtLtI i 


(gcA-^fm-xc 

GGE 


(31) 


,  r  §Cb^3_ 

8ch^2  GGE 


+ 

(32) 


These  equations  may  be  expanded  to  encompass  many  DOC  reservoirs,  but  each  must  be 
defined  by  half-saturation  constant  and  elemental  composition  parameters. 

Nitrate  (+nitrite)  nitrogen  is  budgeted  as  the  difference  between  the  gains  from 
nitrification  ( Nitro )  and  losses  due  to  phytoplankton  uptake: 


8N(X 

dt 


f  3 


Nitro  - 


ZSrjP, 


V  i=i 


(33) 


Nitrification  is  described  by  a  simple  temperature-adjusted  rate: 

Nitro  =  Nre~Kt(T~30)[NH4  ]  (34) 
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The  ammonium  state  equation  is  rendered  more  complex  due  to  residual  nitrogen  fluxes 
from  various  organic  matter  transformations.  The  general  state  equation  is: 


8NH4 

dt 


=  net  _  phytoplankton  +  grazing  _  residuals  +  net  _  bacteria 


+DET  _  residual  +  photochemical  _  gains  -  Nitro 


(35) 


The  first  term  on  the  RHS  is  the  net  phytoplankton  uptake,  which  is  the  difference 
between  phytoplankton  assimilation  and  basal  respiration.  Below  we  show  the  expansion 
of  these  terms  for  Neptune- 1 : 

3 

net  _  phytoplankton  =  V  sx  0[Pi  -  gr  f P0{  ( I  -  Xl  )  (36) 

i=i 


The  second  term  is  the  phytoplankton  grazed  carbon  residuals.  All  of  the  grazed  carbon 
that  is  respired  has  the  appropriate  quantity  of  nitrogen  returned  to  the  ammonium 
reservoir.  The  residual  expression  is  then  the  total  grazed  phytoplankton  carbon  minus 
the  nitrogen  content  transferred  to  the  DOC  and  large/small  detritus  compartments. 

3 

grazing _ residuals  =  V  m]  jP(0l-a2  j0i- a 3  i04-a4  t04-a5  t05)  (37) 

i=i 


The  third  term  is  the  net  heterotrophic  bacterial  component,  the  difference  between 
potential  uptake  and  respiration/grazing  ammonium  gains: 


net _bacteria  =  B02{s2  +m2)-\grbB02  - 


f  grhB  \ 
yGGE  ) 


03\  (38) 


The  detrital  residual,  DET  residual,  term  budgets  the  ammonium  gains  from  the 
solubilization  of  large  and  small  particulate  detritus  to  DOC. 

2 

DET  residual  =  ^  DETki{04  -  XDET03)  (39) 

i=i 


Photochemical  gains  arise  from  the  photolysis  of  both  allochthonous  and  autochthonous 
humic  substances.  The  beta  term  is  the  fraction  of  the  total  photochemical  yield  that  is 
rendered  as  colorless  or  bleached  DOC,  which  is  added  to  the  semilabile  DOC  pool.  The 
remainder  is  remineralized  and  the  nitrogen  content  is  returned  to  the  ammonium 
nitrogen  reservoir. 


2 

photochemical _ gains  =  "Y CDOC nh  i(l-j3i)0s  .  (40) 

i=i 

Here  again,  the  expression  may  be  expanded  to  encompass  many  more  chromophoric 
organic  matter  reservoirs  or  reduced  to  a  single  state  variable. 
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Finally,  nitrification  losses  of  ammonium  to  the  nitrate  compartment,  as  described  by 
equation  (34),  are  deducted  from  the  total.  Alternate  material  pathways,  as  in  Neptune-2, 
for  example,  are  accounted  for  by  revising  the  residual  calculations. 

Dissolved  Inorganic  Phosphorus  (DIP),  largely  as  phosphate,  is  budgeted  analogously  to 
ammonium,  but  without  any  nitrification  terms: 

8DIP  _  _  phytoplankton  +  grazing  _  residuals  +  net  _  bacteria 

ct 

+DET  _  residual  +  photochemical  _  gains  (41) 


For  Neptune- 1,  no  dissolved  organic  phosphorus  (DOP)  is  budgeted:  it  is  assumed  that 
DOP  is  preferentially  remineralized  and  this  process  is  handled  implicitly.  This  renders 
the  photochemical _gains  tenn  in  (41)  as  0  and  the  bacteria  become  competitors  with 
phytoplankton  for  DIP.  However,  by  parameterizing  the  C:P  ratio  of  DOC  to  >  0  the 
model  then  switches  to  potential  net  gains  of  DIP  from  bacterial  organic  matter 
assimilation  and  respiration. 


DOC  for  Neptune- 1  is  produced  via  phytoplankton  grazing,  solubilization  of  particulate 
detritus,  photobleaching  of  CDOC  and  is  lost  via  heterotrophic  bacterial 
assimilation/respiration: 


dDOC 
8 1 


±CDOCpi_A 


+ 


J 

z 


m. 


U,  ,P 


z. 

PetiKt  (1  —  X DET  ) 


g*B  3 

^GGE  J 

(42) 


The  first  two  terms  on  the  RHS  may  be  expanded  or  reduced.  For  Neptune-2,  labile  DOC 
is  produced  exclusively  via  phytoplankton  grazing  and  is  consumed  by  bacteria. 


8DOC \ 
8 1 


3 

Xmi  Ja2_ipi 
i=l 


SrhB 

GGE 


X 


DOC 


(43) 


The  bacterial  uptake  of  DOC  i  not  assimilated  is  transferred  to  both  the  CDOC  and 
semilabile  DOC  pools,  an  additional  parameter,  Xcdoc,  is  required  to  detennine  how  this 
carbon  is  divided  between  the  two  potential  receiving  reservoirs.  Semilabile  DOC  is  thus 
produced  as  a  consequence  of  bacterial  modification  of  labile  DOC,  photobleaching  of 
CDOC,  and  detrital  solubilization.  It  may  also  be  potentially  consumed  by  the  bacteria,  if 
balanced  growth  requirements  are  met. 


dPOC2 

8t 


±CDOClh_A 


+ 


Pet.Kt(\  XDET) 

.  ;=i 

(44) 


+  SrbB(\  ^ CDOC  ) 


GGE 


—  1 
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Small  organic  carbon  particulate  detritus,  SDET,  is  the  sum  of  the  phytoplankton 
production  terms  (via  grazing)  minus  the  solubilization  and  sinking: 


5SDET 

dt 


3 


JP, 


SDETk , 


dwdlSDET 

dz 


(45) 


as  is  large  particulate  detritus,  LDET : 


dLDET 

dt 


3 

ia4_iPi 
i= 1 


LDET  Kt 


dwdlLDET 

dz 


(46) 


The  only  significant  difference  between  how  large  and  small  detritus  are  treated  by  the 
model  is  the  rate  of  sinking.  This  is  because  ecological  models  that  seek  to  mimic  vertical 
carbon  flux  tend  to  assign  high  rates  of  sinking  to  large  detrital  particles  [Walsh,  et  ai, 
1999],  whereas  much  smaller  particles  that  are  still  operationally  defined  as  particulate 
(but  sink  much  more  slowly)  are  important  for  optical  considerations. 

The  detritus  C:N:P  ratio  for  both  large  and  small  detritus  (DET)  is  initially  fixed  at 
250:25:1  reflecting  the  nitrogen  and  phosphorus  depletion  of  non-living  particulate 
organic  matter  in  coastal  systems  relative  to  Redfield  ratios.  Solubilization  and 
remineralization  of  detritus  are  the  two  modeled  loss  pathways,  both  described  by  a  decay 
rate  maximum,  k,„,  at  30°C,  with  the  temperature  dependence  given  by: 

k,  =  min  [Km ,  (/fineA7(r'30) )]  (47) 

Remineralization  of  DET  represents  the  implicit  respiration  of  carbon  by  particle 
adhesive  microbes.  The  biomass  of  particle  adhesive  microbes  is  not  modeled  explicitly, 
but  their  C:N:P  ratios  are  assumed  to  be  similar  to  water-column  heterotrophic 
bacterioplankton— typically  phosphorus  rich  relative  to  carbon  (C:P  ~  53).  The  C:P  ratio 
of  DOC  is  not  applicable  here  since  it  is  assumed  that  DOP  is  preferentially  remineralized 
from  particles,  as  is  often  observed  [Hopkinson,  et  al.,  1997;  Lucea,  et  ah,  2003],  The 
C:P  ratio  of  particle  adhesive  microbes  determines  the  fraction  of  detrital  organic  carbon 
remineralized  to  DIC,  Xdet,  by  calculation  of  the  required  carbon  assimilation  for  particle 
adhesive  microbes: 


DETk& 

4 


DETKtXDET,  DET  =  SDET  +  LDET 


(48) 


The  numerator  on  the  left  hand  side  of  (48)  is  the  total  phosphorus  fluxed  out  of  the  DET 
pool  (detritus  carbon  times  the  rate  of  loss  times  the  P:C  ratio  of  detrital  organic  matter). 
The  assumption  is  made  that  all  of  this  phosphorus  is  used  by  microbes  (the 
remineralization  pathway),  so  if  this  term  is  multiplied  by  a  bacterial  C:P  ratio  (inverse  of 
P:C  for  bacteria,  Eg)  the  total  carbon  flux  to  the  remineralization  pathway  (an  implicit 
overturn  of  the  microbial  population)  is  calculated.  The  expression  may  be  rearranged 
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and  simplified  to  the  ratio  of  the  bacteria  C:P  (53)  to  detritus  C:P  (250)  determining  the 
remineralized  fraction  (53:250  =  21%).  Repeating  the  calculation  for  nitrogen  with  the 
remineralization  fraction  fixed  at  21%  but  allowing  for  a  labile  DOC  pool  C:N  ratio  of  15 
yields  a  particle  adhesive  microbe  C:N  ratio  of  3. 18,  within  the  range  of  observation  for 
marine  bacteria  \Kirchman,  2000], 

Terrestrial  boundary  conditions  serve  as  optional  state  variables.  Terrestrial  CDOC, 
terrestrial  particulate  organic  detritus,  and  suspended  lithogenic  sediments  are  subjected 
to  the  same  advection/diffusion  as  other  state  variables.  The  loss  terms  for  suspended 
sediment  is  sinking  and  the  terrestrial  (or  fluvial)  CDOC  is  subjected  to  photochemical 
decay  which  is  discussed  in  section  2.3.  Organic  detritus  is  subject  to  sinking  and 
solubilization  to  semilabile  DOC. 


Bulk  changes  is  dissolved  inorganic  carbon  (DIC)  are  budgeted,  but  the  model  does  not 
presently  compute  the  carbonate  system  (XDIC,  alkalinity,  pH)  nor  the  potential  air/sea 
gas  exchange.  Loss  of  DIC  is  budgeted  as  phytoplankton  carbon  assimilation  and  gains 
are  a  function  of  the  respiratory,  photochemical,  and  solubilization  terms: 


8DIC 

8t 


3 

Zto -<au  +su 

i= 1 


Sr  i)Pi 


+  B(m2  +  s2 )  +  DET  KtXDET 


grbB( 


1 

GGE 


1) [(Neptune -2  :  (1  - XD0C  j]  + 


ICDOC,,,  , 


0-/0 


+ 

(49) 


The  total  system  carbon,  nitrogen,  and  phosphorus  are  checked  to  verify  the  system 
cycles  elements  conservatively. 

Finally,  for  phytoplankton  functional  groups  where  a  maximum  sinking  speed  is 
designated,  the  realized  sinking  rate  is  a  function  of  the  phytoplankton  biomass.  The 
concept  is  to  simulate  the  formation  of  faster  sinking  cell  aggregates  at  high  biomass 
levels.  This  phenomenon  is  particularly  pertinent  to  diatoms,  or  in  our  size-based 
approach,  the  microphytoplankton.  We  choose  an  expression  that  resembles  the 
Michaelis-Menten  uptake  kinetics  so  that  a  half-saturation  biomass  constant,  nw,  can  be 
designated,  i.e.,  the  biomass  where  the  sinking  rate  is  50%  of  the  maximum  value: 


w 


[p] 


phr  _i 


■  +  tf] 


-  W 


ph 


,(50) 


This  concludes  the  description  and  discussion  of  the  core  ecological  element  cycling 
computations  contained  within  the  module  poseidon.  f90.  There  are  two  fundamental 
quantities  that  must  be  calculated  before  the  above  changes  may  be  computed:  (1)  the 
light-limited  phytoplankton  growth  rate;  and  (2)  the  photochemical  yields  from  CDOM. 
These  quantities  are  computed  in  the  optics  module  apollo.f90  and  are  explained  in  the 
following  sections. 
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2.2.  Bio-Optics  and  Photosynthesis 

The  optical  computations  are  coded  in  the  FORTRAN  90  module  apollo.f90,  which 
contains  the  subroutines  oceanblue  (determines  the  visible  light  attenuation  and  light- 
calculated  phytoplankton  carbon  assimilation  rates),  and  tcdoc  (detennines  ultraviolet 
irradiance  attenuation  and  the  photochemical  yields  from  CDOM).  The  subroutine  angle 
in  the  module  solar.  190  calculates  the  appropriate  irradiance  boundary  condition  and  solar 
zenith  angle  as  a  function  of  position,  day-of-year,  and  time-of-day. 

2.2.1.  Spectral  Absorption  Coefficients 

The  additive  nature  of  water  column  inherent  optical  properties  (IOP’s)  greatly  simplifies 
their  computation  as  a  function  of  carbon  reservoirs.  For  example,  the  total  absorption 
coefficient  may  be  expressed  as  the  sum  of  individual  constituents: 

at(X)  =  aw(X)  +  acdm(k)  +  aph(k)  (51) 

where  the  subscript  w  refers  to  pure  seawater,  cdm  refers  to  colored  detrital  matter,  which 
is  the  sum  of  the  dissolved  and  particulate  non-living  organic  matter  absorption,  and  ph 
refers  to  all  living  phytoplankton  absorption.  All  coefficients  have  units  of  inverse 
meters. 

Absorption  of  light  by  pure  seawater  is  taken  directly  from  published  values  [ Pope  and 
Fry,  1997].  Absorption  by  CDM  is  described  by  adding  the  total  CDOC  and  the 
particulate  detrital  absorptions  together: 


acdmW  = 


Ncdoc 

Nj>ET 

^  aCDOCi  C^) 
i=l 

+ 

^  aDE T_  (A) 

i=i 

(52) 


The  initial  version  requires  the  summation  of  marine  CDOC  and  any  potential  terrestrial 
CDOC  as  well  as  the  total  number  of  particulate  detrital  pools.  For  the  open  ocean 
implementations,  a  minimum  of  1  CDOC  reservoir  and  2  forms  of  particulate  detritus  are 
required. 

The  conversion  of  the  CDOC  reservoirs  to  a  spectral  absorption  value  is  dependent  upon 
multiplication  by  the  mass  specific  absorption,  a*,  parameterized  for  a  reference 
wavelength  (m“  mol'  ;  440  mn).  The  spectral  shape  of  the  mass  specific  absorption  values 
are  determined  by  an  exponential  function  that  is  dependent  upon  a  spectral  slope 
parameter,  S: 


a  *<*  (A)  =  a  (440  nm)  (53) 

For  particulate  detritus,  as  well  as  all  particles  besides  living  phytoplankton  cells,  the 
general  methods  described  by  \Stramski,  et  a!.,  2001],  referred  to  hereafter  as  SBM  01, 
and  the  methods  of  [ Mobley  and  Stramski,  1997],  referred  to  hereafter  as  MS  97,  are 
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Table  4.  Particulate  Absorption  Parameters. 


Component 

Absorption 

units 

ref. 

Organic  Carbon  Detritus 

Ga,da(X)  =  8.791  X  lO'4  exp(-0. 00847  X) 

pm2/particle 

SBM01 

Mineral  Sediment 

ca,min(X)  =  1.013  X  10'3  exp(-0. 00846  X) 

pm2/particle 

SBM01 

Heterotrophic  Bacteria 

°a.hb*c(V  =  -  6.0X10-18(L)  +  5X10-15 

m2/particle 

MS  97 

Detritus  Carbon  g  C  m'3  conversion  to  particles  m'3  =  1.402  x  10'16  g  C 

/  particle 

(Calculated  assuming  a  specific  detritus  scattering  coefficient  of  1.0  m2  g’1  and  an  organic  SPM 
to  organic  carbon  ratio  of  2.6  =  organic  detritus  carbon  specific  scattering  2.6  m2  g  C"1  (555  nm),  then  back 
calculated  using  Stramski  et  al.,  2001  expressions.) 

Mineral  Sediment  g  m'3  conversion  factor  to  particles  m"3  =  3.23  XI O"1^  g  /  particle 

(Calculated  assuming  a  mineral  specific  scattering  coefficient  of  0.5  m2  g'1  (555  nm),  then  back  calculating 

using  Stramski  et  ah,  2001  expressions.) 

Bacteria  moles  nitrogen  m"3  to  bacteria  particles  per  m"3  conversion 

factor  =  2.94  X  10"14  mol  N/particle  (assuming  2.1  fg  N  per  cell  Lee  and  Furhman  1987) 
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Table  5.  Backscattering  Parameters 


Component 

Backscatter 

units 

ref. 

Organic  Carbon  Detritus 

ohhMeta)  =  5.881  X  lO'4  (k)-0'8997 

pm2/particle 

SBM01 

Mineral  Sediment 

ohh.min(l)  =  1.790  XIO'2  (X)-0-9140 

pm2/particle 

SBM01 

Heterotrophic  Bacteria 

Gbb,hbac(K)  =  3.16  x  10-16  for  all  ( X ) 

m2/particle 

MS  97 

Picophytoplankton 

ahhjhytl(X)=U)X  1 0" 1 5  for  all  (X) 

m2/particle 

SBM01 

Nanophytoplankton 

Vbbjhyt2(X)=  1.0  XIO'14  for  all  (X) 

m2/particle 

SBM01 

Microphytoplankton 

VbbjhyaiK)  =  2.0  X  10'13  for  all  (X) 

m2/particle 

SBM01 

Detritus,  Sediment,  and  Bacteria  converted  as  in  Table  4. 

Mass  carbon  m'3  to  particle  m"3  conversion  for  Picophytoplankton 
C:Chl  ratio  =  150.0 

1.433  X  10  03  pg  Chi  per  cell;  Morel  et  al.,  1993 

Mass  carbon  m'3  to  particle  m"3  conversion  for  Nanophytoplankton 
C:Chl  ratio  =  100.0 

1.433  X  10-03  pg  Chi  per  cell;  Bricaud  et  al.,  1988 

Mass  carbon  m'3  to  particle  m"3  conversion  for  Microphytoplankton 
C:Chl  ratio  =  50.0 

1.433  X  10  03  pg  Chi  per  cell;  Ahn  et  al.,  1992 _ 
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utilized  here  to  define  particle  absorption  as  the  product  of  the  particle  concentration,  Np 
(particles  m'  ),  and  the  absorption  cross-section  per  particle,  (7„(a)  (m  per  particle): 


ahbaM)  +  aAM)  +  aseM) 


'NPhbac  ^ a  hbac  W  +  ^dct^a.detWl 
+NP*d°a_rt&) 


(54) 


where  the  subscript  hbac  refers  to  heterotrophic  bacteria,  det  refers  to  particulate  organic 
detritus,  and  sed  refers  to  mineral-based  suspended  sediments.  Table  4  describes  the 
spectral  relationships  for  the  absorption  cross-section  per  particle  and  the  calculated  mass 
concentration  to  particle  density  conversions,  which  are  drawn  from  literature  values  [ Lee 
and Fuhrman,  1987;  Bricaud,  etal.,  1988 ;Ahn,  etal.,  1992;  Morel,  et  al.,  1993], 

The  absorption  coefficient  of  living  phytoplankton  is  calculated  from  the  chlorophyll 
concentration.  Two  steps  required  for  the  computation  are  the  assignment  of  the  carbon  to 
chlorophyll  mass  ratio,  C:Chl,  and  the  spectral  chlorophyll-specific  absorption  coefficient 
a* (a)  (m  mg'  ).  For  the  first  step,  a  static  C:Chl  ratio  is  assigned  to  each  phytoplankton 
functional  group.  For  Neptune- 1  and  Neptune-2,  the  ratios  are  assigned  at  150  for 
picophytoplankton,  100  for  nanophytoplankton,  and  50  for  microphytoplankton. 

Once  the  chlorophyll  reservoir  of  each  phytoplankton  functional  group  is  computed,  it  is 
multiplied  by  the  spectral  chlorophyll-specific  absorption  term  to  compute  the  spectral 
absorption  coefficients.  Variability  exists  within  the  absorption  per  unit  chlorophyll 
parameter  due  to  phytoplankton  taxonomic  differences  related  to  cell  size  and  pigment 
composition  [ Cleveland ,  1995;  Bricaud,  et  al.,  2004;  Millan-Nunez,  et  al.,  2004],  In 
general,  however,  both  field  and  laboratory  culture  observations  confirm  that  a*  often 
decreases  with  increasing  cell  size.  For  the  initial  Neptune  parameterization, 
picophytoplankton  are  given  the  highest  a*  values  and  microphytoplankton  the  lowest. 
The  spectral  values  were  adapted  from  Bricaud  et  al.  [1995]  (see  their  Table  2  and  Figure 
7). 

2.2.2.  Spectral  Backscattering  Coefficients 

Particles  scatter  as  well  as  absorb  incident  irradiance.  Estimating  the  total  backscattering 
coefficient  of  the  water  column  is  essential  to  the  approximation  of  visible  light 
attenuation.  The  backscattering  coefficient  of  pure  seawater  is  assumed  to  be  50%  of  the 
total  scattering  coefficients  across  the  visible  spectrum  (400  -  700  nm),  which  are  taken 
from  literature  values  [Smith  and  Baker,  1981  \  Pope  and  Fry,  1997].  CDOM 
backscattering  is  assumed  to  be  negligible,  leaving  only  particle  backscattering  as 
described  by: 


bbp  (2)  =  bbph  (A)  +  b bhbac  (A)  +  bbdet  (A)  +  bbsed  (A)  (55) 

Each  backscattering  coefficient  may  be  described  as  the  product  of  the  particle 
concentration  and  the  backscattering  particle  cross-section  (m2  per  particle): 
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bb  XX)  = 


NPPh°bb_Ph  W  +  NPhbac°bb  hbac  W  + 


det  W  +  NPSed°bb  sedW 


(56) 


where  ph  refers  to  phytoplankton,  and  the  remaining  subscripts  refer  to  heterotrophic 
bacteria,  organic  carbon  detritus,  and  mineral-based  suspended  sediments,  respectively. 

The  expressions  describing  the  value  of  each  backscattering  cross-section  as  a  function  of 
wavelength  (nm)  is  proved  in  Table  5,  as  well  as  mass  concentration  to  particle  density 
conversion  factors  for  the  three  classes  of  phytoplankton  in  Neptune- 1  and  Neptune -2. 

2.2.3.  Irradiant  Boundary  Condition 

Incident  irradiance  is  modeled  using  the  Gregg  and  Carder  clear  sky  irradiance  model 
[ Gregg  and  Carder,  1990]  for  the  visible  range  (400  -  700  nm),  and  equations  adapted 
from  the  GCSOLAR  UV  model  [. Zepp  and  Cline,  1977]  for  the  lower  spectral  range  (290 
-  400  nm).  The  models  were  combined  into  a  single  source  code  that  gives  direct  and 
diffuse  clear  sky/band-centered  irradiances  just  below  the  sea  surface  with  AX  =  10  nm 
for  band-centered  wavelengths  295  -  555  nm  and  AX  =  20  nm  for  band-centered 
wavelengths  570  to  590  nm. 

A  spectrally  neutral  reduction  factor  of  15%  is  applied  to  the  direct  and  diffuse  output  to 
approximate  the  average  reduction  by  clouds.  Previous  versions  of  this  model  [Jolliff, 
2004]  used  National  Center  for  Environmental  Prediction’s  (NCEP)  reanalysis 
climatology  cloud  fractions  to  better  estimate  the  impact  of  clouds.  This  technique  can  be 
implemented  for  any  specific  area  where  NCEP  estimates  are  available. 

2.2.4.  Attenuation  and  the  Visible  Photon  Budget 

The  incident,  planar,  downwelling,  and  band-centered  irradiance  just  below  the  air/sea 

2  1 

interface,  Ed(0",X;  W  m'  nm'  )  is  first  converted  to  the  band-averaged  photon  flux  (p,  mol 
2  1 

photons  m'  s'  ): 


p,„,W  =  Ed  Jlra) 


AA("m)  (57) 
q{X)a 


and 


P diff  (X)  ~  Eddiff{X) 


AA(nm) 

q(X)a 


(58) 


where 


q(A)  = 


he 

A(m) 


(59) 
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and  h  is  Planck’s  constant,  c  is  the  speed  of  light,  and  X  is  wavelength  in  meters.  The  a  in 
the  RHS  of  equations  (57)  and  (58)  is  Avogadro’s  number  and  the  wavelength  bandwidth 
in  these  equations  is  in  nanometers. 

The  photon  flux  is  then  progressively  attenuated  over  descending  depth  increments,  Az, 
by  estimation  of  the  downwelling  attenuation  coefficients  for  the  direct  and  diffuse 
photon  fluxes,  Kd  dufA)  and  Kd  dm(A),  using  the  single-scattering  approximation  [ Gordon , 
1989]: 


Kd  dir(z^) 


aXz,X)  +  bb,(z,X) 

'  0, 


and 


at(z,A)  +  bbt(z,X) 

0.7  1  ' 


where  the  direct  attenuation  is  divided  by  the  estimated  direct  path  length,  i.e.,  the  cosine 
of  the  solar  zenith  angle,  0Z,  once  refracted  through  the  air/sea  interface.  The  diffuse 
attenuation  is  normalized  by  the  cosine  of  a  constant  diffuse  angle  of  -45°. 

Attenuation  over  each  discrete  vertical  depth  interval  is  calculated  following  the  Beer- 
Lambert  Law: 


Pd  dir(Z’  Pd  dir 


(z-  \,A)e[~Kd-di' 


( z,A)dz ] 


(62) 


and 


(63) 


The  total  quantity  of  photons  absorbed  by  the  phytoplankton  absorption  fraction  is 
retained  over  each  depth  level  and  used  to  calculate  the  light-limited  growth  rate  and 
photoinhibition  terms.  Since  backscattered  photons  may  be  potentially  absorbed  from  the 
upwelling  irradiance  stream,  some  estimation  of  this  contribution  is  required  to  better 
estimate  the  total  absorbed  photon  budget.  Photons  absorbed  from  the  upwelling 
irradiance  stream  are  thus  estimated  by  calculating  the  irradiance  reflectance  (R  =  Eu  /  Ed) 
for  the  discrete  depth  level  [ Mobley ,  1994]: 


R(z,X) 


(A/z,T)-q,(z,A)///1) 
(Ku  (z,  2)  +  cit  (z,  X)/i~v) 


where  the  average  cosine  of  upwelling  irradiance,  //u,  is  held  constant  at  -0.39  [Kirk, 
1994].  The  diffuse  and  direct  notation  was  dropped  for  brevity;  for  the  direct  stream,  /id  = 
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0Z,  and  for  the  diffuse  stream  p d  =  0.7.  Since  the  attenuation  functions  are  approximated 
using  the  single  scattering  approximations,  equation  (64)  may  be  recast  as: 

R( z,X)  =  r - ,  ,  (65) 

[2  at(z,A)  +bbt(z,A)\jUu 

The  upwelling  stream  for  each  depth  interval  is  the  downwelling  flux  multiplied  by  the 
irradiance  reflectance. 

The  important  quantity  towards  which  we  are  calculating  is  the  photon  flux  absorbed  by 
phytoplankton  over  each  discrete  depth  increment,  Az.  We  begin  by  first  tallying  the 
downwelling  photon  flux  density  absorbed  by  phytoplankton,  pdv( z,  X),  then  calculating 
the  upwelling  analogue  quantity,  puv( z,  X): 

Pdv  ( Z >  V  =  Q([/7/_^  (Z  ~hV~  Pd_dir  (Z>  V]  +  [ Pd_diff  (Z  -IV-  Pd_diff  V’  'O])  (66) 


and 


Q  = - - -  (67) 

at  (z,  /t)  +  bbt  (z,  A)  Az 

Upwelling  photons  are  also  attenuated,  although  the  upwelling  irradiance  increases  with 
ascending  depth.  We  thus  approximate  a  discrete  depth  upwelling  attenuation  solely  as  a 
function  of  the  depth-centered  irradiance  reflectance: 

pm(z,A)  =  n(RJlf/(z,Mpi  m(z,X)  +  Rllr(z,X)pd  *(^>)(1-e'"w|1'1)(68) 

The  spectrally  integrated  total  photon  flux  harvest  for  phytoplankton,  nH  (mol  photons 
m‘  s"  )  between  grid  levels  z-1  and  z  is  then  the  spectral  summation  of  the  calculations 
centered  on  depth  level  z: 


N 

pH  =  Yj  Puv  V,  M )  +Pdv  V,  A  A, )  (69) 

i= 1 

From  first  principles,  only  photons  absorbed  are  available  to  do  photochemical  work.  For 
the  enzymatically  mediated  photochemical  process  of  photosynthesis,  this  simple 
relationship  can  be  expressed  through  the  apparent  quantum  yield: 

d>c  =  mole  organic  carbon  produced  /  mole  photons  absorbed  (70) 

The  theoretical  phytoplankton  carbon  yield  (TCPY;  mol  carbon  in’  s’  )  is  then  simply  the 
calculated  photon  harvest  multiplied  by  the  apparent  photosynthetic  quantum  yield: 
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TPCY  =  pH  (®c)  (71) 


The  maximum  rate  of  photon  absorption  for  the  phytoplankton  biomass  (pM;  mol 
photons  m"  s'  )  may  then  also  be  calculated  as  the  product  of  the  maximal  temperature- 
adjusted  rate  of  growth,  the  phytoplankton  biomass,  and  the  inverse  of  the  quantum  yield: 

pM  =  g^P{®cy  (72) 

It  follows  that  absorbed  photon  fluxes  below  this  threshold  would  limit  growth,  while 
absorbed  photon  fluxes  above  this  threshold  may  potentially  inhibit  growth.  However,  a 
strictly  linear  percentage  reduction  in  maximal  growth  rate  does  not  reproduce  the  normal 
hyperbolic  relationship  observed  between  growth  (or  photosynthetic  rate)  and  irradiance 
(i.e.,  the  PI  curve).  This  is  because  photosynthesis  is  an  enzymatically-mediated  process, 
and  the  kinetics  of  such  processes  is  best  mimicked  by  Michealis-Menten-type  functions. 

If  the  ratio  of  TPCY  to  the  temperature-adjusted  maximal  assimilation  rate  (gmaxi)  times 
the  phytoplankton  biomass  (P)  is  less  than  unity,  the  light-limited  growth  calculation  is 
performed: 


Sl  = 


¥ 

MM  p  +  ¥ 


(g 


maxT 


) 


(73) 


where 


¥ 


TPCY 

&  max  T  ^ 


(74) 


and  MMP  is  the  Michaelis-Menten  half- saturation  constant  for  photon  flux,  empirically 
detennined  to  be  ~  0.072  for  a  wide  variety  of  cultured  phytoplankton  [Rubio,  et  al., 
2003]. 

If  the  TCPY  ratio,  ¥,  is  greater  than  4.0  then  the  estimated  photoinhibitive  effect  on 
growth  rate  is  estimated  as  a  function  of  the  incident  photon  flux  [Rubio,  et  al.,  2003]: 


Sl  = 


1.0 


1.0  + 


20-0  yjpdyz-\) 


( g max  T  ) 


(75) 


and 

PdAZ-l)  =  Y,[pd  dir(Z-1’4)  +  Pd  ^(z-1,4)]  (76) 

i= 1 
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These  expressions  mimic  the  photosynthesis-irradiance  (PI)  curve  of  marine 
phytoplankton,  but  are  instead  calculated  as  a  function  of  the  phytoplankton  community’s 
maximal  rate  of  growth,  biomass,  quantum  yield  of  photosynthesis,  and  the  absorbed 
photon  flux. 

2.3.  Ultraviolet  Attenuation  and  Photochemical  Sub-Model. 

The  attenuation  of  ultraviolet  light  and  the  computation  of  photochemical  yields  follow  a 
construct  similar  to  the  one  outlined  in  Section  2.2:  the  attenuation  of  photons  with  depth 
is  budgeted  such  that  the  photon  quantity  absorbed  by  the  reactive  constituent  is  retained 
and  multiplied  by  a  quantum  yield  to  determine  the  product  yield  over  each  discrete  depth 
increment.  The  boundary  visible  irradiance  values  were  calculated  as  a  function  of  the 
solar  zenith  angle  which  was,  in  turn,  calculated  as  a  function  of  latitude,  longitude,  day- 
of-year,  and  time-of-day  using  equations  set  forth  in  Gegg  and  Carder  [1990],  The 
expressions  that  describe  incident  downwelling  ultraviolet  irradiance  as  a  function  of  the 
solar  zenith  angle  were  adapted  from  the  GCSOLAR  model  version  1 .2 
(http://www.epa.gov/ceampubl/swater/gcsolar/)  wherein  the  computations  are  based  on 
the  work  of  Baker  et  al.  [1980],  The  UV  spectra  are  centered  on  the  295  -  395  nm  bands 
at  Aa=  1 0  nm  resolution. 

The  stratospheric  ozone  amount  must  be  designated  before  determining  the  UV  boundary 
condition.  Either  a  default  value  may  be  used  or  the  UV  boundary  condition  may  be 
determined  from  Total  Ozone  Mapping  Spectrometer  global  data  (courtesy  of 
NASA/GSFC;  http://toms.gsfc.nasa.gov/ozone/ozoneother.html).  For  example,  the  grand 
mean  of  the  zonal  averages  from  each  latitude  band  from  the  entire  time  series  (1997  - 
2005)  was  used  determine  the  ozone  value  for  the  Gulf  of  Mexico  and  the  Sargasso  Sea. 

The  first  1 1  spectral  bandwidths  of  the  irradiance  boundary  encompass  the  ultraviolet 
range  290  -  400  nm.  The  diffuse  and  direct  components  of  the  planar,  downwelling,  and 
band-centered  irradiance  just  below  the  sea  surface,  Ed(0‘,X;  W  m'  nm'  )  are  first 
converted  to  the  band-averaged  photon  flux  (p ,  mol  photons  m'  s'  )  following  equations 
(57-59). 

We  make  the  simplifying  assumption  that  scattering  within  the  ultraviolet  (UV)  range  is 
negligible  compared  to  absorption,  such  that  the  ultraviolet  attenuation  functions  simplify 
to: 


Kt  =  (77) 

O- 


and 


■)=  a'(Z0' j“,'>  (78) 
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where  the  direct  attenuation  is  divided  by  the  estimated  direct  path  length,  i.e.,  the  cosine 
of  the  solar  zenith  angle,  0Z,  once  refracted  through  the  air/sea  interface.  The  diffuse 
attenuation  is  normalized  by  the  cosine  of  a  constant  diffuse  angle  of  -45°. 

We  further  assume  that  CDM  and  water  dominate  the  UV  absorption: 

Ct-tty uv)  =  Clw(k\jv)  +  ctedmO^uv)  (79) 

Absorption  of  light  by  pure  seawater  is  taken  directly  from  published  values  and 
absorption  by  CDM  is  described  by  adding  the  total  CDOC  and  the  particulate  detrital 
absorptions  together: 


cdm 


IV  IV 

aCDOC ,  (A/F  )  aDET,  (■ A/F  ) 


(80) 


Attenuation  of  the  downwelling  UV  photon  flux  is  calculated  as  in  equations  (62-63).  No 
reflectance  is  considered  since  scattering  is  assumed  to  be  negligible. 


The  UV  photon  harvest  for  each  discrete  depth  increment  by  CDOC  component  being 
considered  is  then  described  by: 


PdV(Z’  A/f)  =  q 


\_Pd_dir(Z  1’A Jv)  Pd_diriZ’  A/f)]  + 
\Pd_diff  (Z  _  f  A JV  )  _  Pd  diff  ( Z ’  A/F  )] 


(81) 


and 


Q 


aCDM\Z’>\jy)  ^ 

®t(Z’  A/f)  ^ z 


(82) 


The  quantum  yields  for  photochemical  products  from  CDOM  are  adapted  from 
Johannessen  and  Miller  [2001]  for  the  CO?  yield.  They  give  a  coastal  value: 

CO2(A)  =  e-(636+0-0140(i-290))(83) 


and  an  open  ocean  value: 


Oph  ca(T)  =  U(5'53+000914U-290,)(84) 

which  we  apply  to  terrestrial  CDOC  and  marine  CDOC,  respectively.  We  assume  that  the 
yield  of  CO?  from  CDOM  photochemical  change  is  roughly  equal  to  the  yield  of 
bleached  DOC  [ Miller  and  Moran,  1997;  Mopper  and  Kieber,  2000],  such  that  the  total 
loss  of  CDOM  from  photochemical  change  is  double  the  spectrally  integrated  product 
yield  of  CO?: 
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11 

total _ph_yield(z)  =  2.0 YjPdV(zAvi)°ph_co2  ('UX85) 

1=1 

Thus  half  of  the  total  photochemical  yield  is  diverted  to  the  colorless,  semilabile  DOC 
pool  and  half  returns  to  the  DIC  pool. 

Previous  application  of  the  photolysis  sub-model  [. Jolliff,  et  al.,  2003;  Jolliff,  2004] 
suggested  that  CDOC  likely  consists  of  a  photochemically  reactive  component  and  an 
inert,  or  background  component.  We  thus  modify  the  total  absorption  description  to 
accommodate  a  background  absorption  CDM  signal,  acDMb,  that  is  not  linked  to  any  state 
variable  carbon  reservoir: 


tb(Xjv)  —  ClwQ\. uv)  T  O-cdmQ^ Uv)  T  tlcdmbO^Uv}  (86) 


where  the  background  signal  is  a  constant  0.015  nf1  at  440  nm,  and  is  spectrally 
described  by  equation  (53)  using  a  spectral  slope,  S,  of  0.015.  Application  of  this 
formulation  to  the  Sargasso  Sea  suggests  the  background  signal  is  closer  to  0.006  m'1  at 
440  nm.  Discerning  ways  to  represent  photochemically  active  and  inert  chromophores  in 
the  model  remains  an  active  area  of  research  and  development. 

2.4.  MOD AS-Derived  Vertical  Diffusivity  for  the  One-Dimensional  Case 

For  the  one-dimensional  application  of  the  model,  the  eddy  mixing  coefficients,  Kz(z), 
may  be  inferred  from  the  Modular  Ocean  Data  Assimilation  System  vertical  temperature 
fields  without  any  additional  information.  To  accomplish  this  task,  we  adapt  the  vertical 
mixing  parameterization  scheme  of  Pacanowski  and  Philanderer  [1981],  or  PP81  for 
brevity.  Their  scheme  solves  for  Kz(z)  as  a  function  of  the  local  gradient  Richardson 
number: 


Kz{z )  = 


(i+«Rigr+i 


+  Kb  (87) 


2  1  2  1 
where  vo=  100  cm  s'  ,  Kb,  the  background  diffusivity,  is  0.1  cm"  s'  ,  and  the  empirical 

constants  are  assigned  the  PP81  values  of  n  =  2,  and  a  =  5.  As  in  Li  et  al.  [2001],  the 

local  Richardson  gradient,  Rig,  is  calculated  as  the  ratio  of  the  Brunt-Vaisala  frequency  to 

the  vertical  current  shear.  The  Brunt-Vaisala  frequency  is  calculated  from  the 

temperature-based  density  gradients  (assuming  a  constant  salinity)  from  solution  of  the 

UNESCO  equation  of  state.  Since  the  vertical  current  shear  is  unknown,  we  back- 

calculated  a  bulk  shear  approximation  by  assuming  the  vertical  eddy  diffusivity  at  depths 

above  the  isothennal  layer  depth  are  —100  cm2  s'1.  The  isothermal  layer  depth  is  defined 

as  the  first  descending  depth  below  10-meters  where  the  temperature  difference  from  10- 

meters  depth  exceeds  0.8  deg  C  [Kara,  et  al.,  2000],  The  mean  Brunt-Vaisala  within  the 

isothermal  layer  depth  divided  by  the  maximal  diffusivity  yielded  a  bulk  shear  term  that 

was  used  to  calculate  eddy  diffusivities  below  the  isothermal  layer  depth  using  the  PP8 1 

scheme. 
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2.5.  Model  Integration 

The  flow  of  the  code  execution  is  shown  in  Figure  4.  The  main  program  unit 
declares/initializes  the  state  variables,  and  then  opens/reads  boundary  conditions  and 
hydrodynamic/temperature  files  necessary  for  code  execution.  The  irradiant  boundary 
condition  is  computed  ahead  of  code  execution  for  a  range  of  solar  zenith  angles,  and 
then  selected  during  code  execution  depending  upon  position  and  time.  The 
photochemical  yield  calculations  (Section  2.3)  and  light-limited  growth  calculations 
(Section  2.2)  are  computed  first  and  the  necessary  information  (light-limited  growth  rate, 
photochemical  yields)  is  passed  to  the  element  cycling  subroutine  (Section  2.1).  The 
hydrodynamic  advection/diffusion  and  particle  sinking  computations  are  performed  and 
then  the  code  loops  to  another  time-step  execution.  The  state  variables  as  well  as  some 
optical  infonnation  are  recorded  at  the  end  of  each  daily  iteration. 

For  the  one-dimensional  application,  the  MOD  AS  script  to  generate  temperature  fields 
and  the  FORTRAN  code  to  generate  eddy  diffusion  coefficients  must  be  run  ahead  of 
model  execution.  Presently,  we  are  combining  5-years  or  more  of  daily  MOD  AS  fields 
into  a  composite  climatology  for  a  given  area  of  interest.  The  MODAS  script,  the  eddy 
diffusion  and  temperature  climatology,  and  the  main  program  unit  must  all  generate  the 
same  grid  dimensions  and  resolution.  The  latitude  and  longitude  location  must  also  be 
specified  in  the  main  program  unit  in  order  to  correctly  calculate  the  sun  angle  as  a 
function  of  time  and  location. 
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Figure  4.  Model  code  execution  flow  chart. 


Jolliff  and  Kindle 


Page  34  of  49 


END 


Figure  4.  Continued,  Model  code  execution  flow  chart. 
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3.  Example  Application:  Western  Gulf  of  Mexico 

Briefly,  we  provide  an  example  of  our  application  of  Neptune  to  the  western  Gulf  of 
Mexico.  The  MOD  AS  script  was  used  to  generate  daily  temperature  fields  for  five  years 
from  1  to  161-meters  depth  at  1 -meter  resolution.  These  data  were  merged  into  a  single 
annual  climatology  (Figure  5)  and  the  eddy  diffusion  coefficients  were  estimated  from 
this  field.  The  nutrient  fields  (NO3,  NH4,  DIP)  were  initialized  using  temperature  to 
nutrient  relationships  observed  in  the  northeastern  Gulf  of  Mexico  [. Jochens ,  et  al.,  2002]. 
The  model  was  then  run  towards  a  steady  state  solution:  the  annual  cycle  was  repeated 
for  five  years  and  the  system  was  forced  to  element  conservation  by  implicitly 
remineralizing  particulate  detritus  within  the  lower-most  grid  cell. 

Our  detailed  model  description  of  the  relationship  between  carbon  reservoirs  and  optical 
properties  allowed  us  to  compare  the  steady  state  model  solution  with  Sea-viewing  Wide- 
Field-of-view  Sensor  (SeaWiFS)  mean  climatology  extracted  from  a  1°  x  1°  area 
surrounding  our  location  of  interest.  Both  the  standard  OC4v4  chlorophyll  algorithm  and 
the  QAA  algorithms  [Lee,  et  al.,  2002]  were  applied  to  the  extracted  SeaWiFS  data  so 
that  model/satellite  comparisons  of  surface  phytoplankton  absorption  (443  mn),  surface 
chlorophyll  concentration,  and  CDM  absorption  (412  nm)  could  be  made  (Figure  6).  The 
model  mimics  both  the  amplitude  and  phase  of  the  satellite  signals.  We  were  able  to  use 
the  satellite  data  to  begin  to  constrain  model  structure  and  model  parameter  selection,  as 
described  elsewhere  [. Jolliff,  et  al.,  2006].  The  model  also  reproduced  the  depth 
distribution  of  chlorophyll  common  to  stratified  waters,  i.e.,  the  appearance  of  the  deep 
chlorophyll  maximum,  particularly  during  summer  months  (Figure  7). 

An  interesting  result  of  our  model  is  that  the  DOC  and  CDOC  distributions  are  out  of 
phase  (Figure  8).  In  fact,  DOC  apparently  accumulates  in  surface  waters  during  summer 
whereas  CDOC  is  depleted  during  this  time.  This  type  of  CDOC  to  DOC  decoupling  has 
been  observed  in  subtropical  oceans  [Nelson,  et  al.,  1998;  Steinberg,  et  al.,  2004].  Our 
continuing  sensitivity  analyses  further  demonstrate  that  the  cycling  of  DOC  impacts  the 
recycling  of  nutrients  and,  in  turn,  the  amplitude  and  phase  of  the  seasonal  surface 
chlorophyll  signal.  We  are  now  investigating  these  model  features  using  in  situ  data  from 
the  Sargasso  Sea.  The  portability  of  the  combined  Neptune/MODAS  modeling  construct 
allows  us  to  examine  how  data  may  constrain  the  model  in  any  region  around  the  globe 
where  data  are  available.  We  also  plan  to  couple  refined  versions  of  Neptune  to  regional 
applications  of  the  Naval  Research  Faboratory  Coastal  Ocean  Model  so  as  to  apply  our 
methods  of  satellite  validation  and  constraints  to  three-dimensional  fields. 
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Figure  5.  (A)  Map  of  western  Gulf  of  Mexico  test  area  for  Neptune,  and  (B)  the  5-year 
annual  MODAS  climatology  for  one-dimensional  temperature  profiles. 
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Figure  6.  The  model’s  surface  results  are  compared  to  the  mean  annual  cycle  (2  years  are  shown)  as  estimated  from  daily  averages  of 
SeaWiFS  data  (2002  -  2004)  of  phytoplankton  absorption  (A),  surface  chlorophyll  concentration  (B),  and  CDM  absorption  (C). 

The  satellite  data  suggest  that  the  relationship  between  chlorophyll-a  and  phytoplankton  absorption  is  not  constant;  phytoplankton 
absorption  efficiency  per  unit  chlorophyll-a  (a*  m2  mg-1  pigment)  decreases  during  the  winter.  The  model  mimics  this  pattern  via 
changes  in  the  phytoplankton  composition:  during  winter,  microphytoplankton  are  more  abundant.  Microphytoplankton  are 
parameterized  with  a  lower  a*  than  picophytoplankton. 


Figure  7.  The  model  simulates  the  transitions  between  the  winter/spring  well-mixed  chl-a  distribution  and  the  summer  stratified  water 
column  wherein  a  deep  chlorophyll  maximum  is  evident.  Simulated  surface  waters  during  summer  are  dominated  by  the 
picophytoplankton  size  class. 


DOC  (no  color)  ( xM  C  Colored  DOC  | xM  C  Heterotrophic  Bacterial  Biomass  \xM  C 


Figure  8.  The  model  also  shows  opposite  trends  for  colorless,  bulk  dissolved  organic  carbon  (DOC)  and  colored  dissolved  organic  carbon 
(CDOC)  distributions:  DOC  accumulates  in  surface  waters  during  summer  whereas  surface  water  CDOC  is  depleted  during  this 
time.  Accumulated  DOC  is  eventually  broken  down  by  a  heterotrophic  bacterial  biomass  that  peaks  during  late  summer,  early  fall. 
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4.  Summary 

A  modeling  system  has  been  constructed  that  combines  ecological  element  cycling, 
photochemical  processes,  and  bio-optical  processes  into  a  single  simulation  that  may  be 
coupled  to  hydrodynamic  models  that  provide  temperature  fields  as  well  as  the 
advection/diffusion  of  state  variables.  The  model  is  derived  from  a  history  of  ocean 
biogeochemical  models  that  describe  the  transfonnation  of  elemental  reservoirs  (carbon, 
nitrogen,  phosphorus)  based  upon  lower  trophic-order  ecosystem  function.  Here  we  relate 
these  reservoirs  to  optical  properties  in  order  to  describe  how  the  spectrally  decomposed 
attenuation  of  both  visible  and  ultraviolet  irradiance  stimulates  photochemical  and 
photosynthetic  changes  in  upper  ocean  elemental  reservoirs.  For  our  initial  model 
startup,  we  have  parameterized  phytoplankton  functional  groups  based  on  size  since  size 
appears  to  impact  both  nutrient  affinity  and  optical  property  variability  in  consistent 
ways. 

Our  model  description  of  the  relationship  between  reservoirs  of  organic  carbon  and 
optical  properties  allows  us  to  use  satellite  ocean  color  data  to  begin  to  validate  and 
constrain  the  model.  We  have  thus  far  coupled  the  model  to  the  Modular  Ocean  Data 
Assimilation  System,  which  permits  us  to  examine  the  one-dimensional  case  in  any 
region  of  interest  around  the  globe.  This  portability  will  further  allow  us  to  continue  to 
refine  and  constrain  the  model  using  satellite  ocean  color  data  as  well  as  in  situ  data, 
wherever  it  becomes  available. 
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Appendix.  Comprehensive  List  of  Symbols  and  Acronyms 
Table  Al.  Ecological  Symbols  and  Acronyms 


Ks 

n 

tiw 

fJ-net 

gr 

gL 

gN 

gr 

.gmaxT 
,§m30 
A netb 
grb 

gCb 

gNb 

gPb 

gmaxTb 

gm30b 

51 

5 2 
Sim 
S2m 
mi 
m2 

m  1  max 
m-2max 
m  1  maxT 
nJ2maxT 

b 

Kt 

Km 

X 

Xdoc 

Xdet 

XcdOC 

MBNA 

BNR 

MBPA 

BPR 

DOC 

CDOC 

T 

Kt 


general  form  of  the  Michaelis-Menten  half-saturation  constant 

model  Michaelis-Menten  (MM)  half-saturation  constants  (see  parameters  list) 

3 

half-saturation  constant  for  phytoplankton  sinking;  mmol  C  m" 
net  community  carbon  assimilation  rate  for  phytoplankton;  s'1 
maximal  realized  gross  carbon  ass’n  rate;  s'1 
light-limited  maximal  gross  carbon  ass’n  rate;  s'1 
nitrogen-limited  maximal  gross  carbon  ass’n  rate;  s'1 
phosphorus-limited  maximal  gross  carbon  ass’n  rate;  s'1 
temperature  adjusted  gross  maximal  carbon  ass’n  rate;  s'1 
maximum  gross  carbon  ass’n  rate  at  30°  C;  s'1 
net  community  carbon  ass’n  rate  for  heterotrophic  bacteria;  s'1 
maximal  realized  net  cellular  carbon  ass’n  het.  bacteria;  s'1 
maximal  carbon-limited  bacteria  net  cellular  carbon  ass’n;  s'1 
maximal  nitrogen-limited  bacteria  net  cellular  carbon  ass’n;  s'1 
maximal  phosphorus-limited  bacteria  net  cellular  carbon  ass’n;  s'1 
temperature-adjusted  gross  cellular  maximal  carbon  assimilation  rate;  s'1 
maximum  gross  cellular  carbon  assimilation  rate;  s'1 
temperature-adjusted  phytoplankton  basal  respiration  rate;  s'1 
temperature-adjusted  bacterial  basal  respiration  rate;  s'1 
maximum  phytoplankton  basal  respiration  rate  at  30°  C;  s'1 
maximum  bacterial  basal  respiration  rate  at  30°  C;  s'1 
realized  phytoplankton  mortality  rate  (grazing);  s'1 
realized  heterotrophic  bacteria  mortality  rate  (grazing);  s'1 
phytoplankton  mortality  rate  (grazing)  at  30°  C;  s'1 
heterotrophic  bacteria  mortality  rate  (grazing)  at  30°  C;  s'1 
temperature-adjusted  phytoplankton  mortality  rate  (grazing);  s'1 
temperature-adjusted  heterotrophic  bacteria  mortality  rate  (grazing);  s'1 
Nitrate  uptake  inhibition  term 

Temperature-adjusted  particulate  organic  detritus  solubilization  rate;  s'1 

Maximum  particulate  organic  detritus  solubilization  rate  at  30°  C;  s'1 

Nitrate  to  total  nitrogen  uptake  fraction 

Fraction  of  Bacterial  DOC  uptake  from  Labile  DOC 

Fraction  of  solubilized  DET  remineralized 

Fraction  of  labile  carbon  transferred  to  CDOC  pool  (Neptune -2) 

Maximum  potential  Bacterial  Nitrogen  Assimilation;  mmol  N  m'3 

3 

Bacterial  Nitrogen  Requirement;  mmol  N  m' 

Maximum  potential  Bacterial  Phosphorus  Assimilation;  mmol  P  m'3 

3 

Bacterial  Phosphorus  Requirement;  mmol  P  m" 

Dissolved  Organic  Carbon;  mmol  C  m'3 

3 

Colored  Dissolved  Organic  Carbon;  mmol  C  m" 

Temperature; 0  C 

Temperature  inhibition  coefficient 
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Nr 

Iv 

a 

e 

£ 

GGE 

B 

P 

N03 

nh4 

DIP 

Nitro 

DET 

LDET 

SDET 

wdi 

Wd2 

Wd3 

Wph 

Wphr 

P 

Vm 

t 


Table  A2. 


X 

a 

h 

c 

pO-) 

Pd  dir(X) 
PddifjQ 0 
Pdv(X) 

Puv(X) 

pH 

pM 

PdX 

Ed  (X) 
Eu  ( X ) 

q 

bb(X) 


Maximum  nitrification  rate  at  30°  C;  s'1 

3  1 

Ivlev  grazing  attenuation  parameter;  m  (mmol  C)' 

Grazing  transformation  fractions  (see  parameters  list) 

Molar  N:C  ratio  (see  parameters  list) 

Molar  N:P  ratio  (see  parameters  list) 

Heterotrophic  Bacteria  Gross  Carbon  Growth  Efficiency 
Heterotrophic  Bacteria  Carbon;  mmol  C  m'3 
Phytoplankton  Carbon;  mmol  C  m'3 

3 

Nitrogen  as  nitrate  (+nitrite);  mmol  N  m' 

3 

Nitrogen  as  ammonium;  mmol  N  m" 

Dissolved  Inorganic  Phosphate  (largely  as  orthophosphate);  mmol  P  m'3 
Quantity  of  NH4  to  N03  conversion  per  unit  time;  mmol  N  m'  s' 
particulate  organic  detrital  carbon;  the  sum  of  LDET  and  SDET ;  mmol  C  m"3 

3 

large  particulate  organic  detrital  carbon;  mmol  C  m" 
small  particulate  organic  detrital  carbon;  mmol  C  m"3 
Sinking  velocity  small  particulate  organic  detritus;  m  s'1 
Sinking  velocity  large  particulate  organic  detritus;  m  s'1 
Sinking  velocity  of  suspended  mineral  sediments;  m  s'1 
Maximum  sinking  velocity  of  phytoplankton;  m  s'1 
Realized  maximum  sinking  velocity  of  phytoplankton;  m  s'1 
Fraction  of  photochemical  CDOC  yield  as  bleached  DOC 
Maximum  uptake  rate  for  ideal  Michaelis-Menten  kintetics 
time;  seconds 


Optical  Symbols  and  Acronyms 

wavelength;  nanometers  unless  otherwise  specified 

Avogadro’s  number 

Planck’s  Constant 

speed  of  light;  m  s'1 

photon  flux;  mol  photon  in'"  s' 

2  1 

downwelling,  direct  photon  flux;  mol  photon  m'  s' 

2  1 

downwelling,  diffuse  photon  flux;  mol  photon  m‘  s' 

3  1 

downwelling  photon  flux  density;  mol  photon  m'  s' 

upwelling  photon  flux  density;  mol  photon  in'  s' 

spectrally  integrated  photon  flux  density  absorbed  (harvested)  by 

phytoplankton;  mol  photon  m'  s' 

maximum  potential  (saturating)  phytoplankton  absorbed 

photon  flux  density;  mol  photon  m'  s' 

spectrally  integrated  sum  of  direct  and  diffuse  incident  photon  fluxes; 
mol  photon  m'2  s'1 

band-centered  downwelling  planar  irradiance  W  m'  nm~ 

2  1 

band-centered  upwelling  planar  irradiance  W  m'  nm~ 

Energy;  Joules 

spectral  backscattering  coefficient;  m'1 
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a(k) 
a  *(A) 
a*  chi(^) 

<2*cdm(k) 

o 

Oc 

®ph 

®ph  C02  0-) 

0z 

Np 

Obb(k) 

TCPY 


'P 

Pd 

Pa 

MMP 
Kd_dir(X) 
Kd_diff  (X) 
Ku(k) 

R(X) 

z 

Q 


spectral  absorption  coefficient;  m'1 

spectral  mass  specific  absorption;  m'1  (mass  units)'1 

chlorophyll-specific  absorption;  m  mg' 

2  1 

cdm  mass  specific  absorption;  m“  mol" 

apparent  quantum  yield;  mol  product  yield  per  mol  photon  absorbed 
apparent  quantum  yield  for  photosynthesis;  mol  carbon  assimilated  per  mol 
photons  absorbed 

apparent  quantum  yield  for  CDOC  photochemistry;  mol  CDOC  altered  per 
mol  photons  absorbed 

spectral  apparent  quantum  yield  for  CDOC  photochemistry  specific  to  CO2 
yield;  mol  CO2  per  mol  photons  absorbed 
Cosine  of  the  refracted  solar  zenith  angle 

■j 

particle  concentration;  particles  m" 
backscattering  particle  cross-section  (m“  particle"  ) 
absorption  particle  cross-section  (m  particle"  ) 

Theoretical  Photosynthetic  Carbon  Yield,  the  product  of  the  absorbed  photon 
flux  density  and  the  apparent  quantum  yield  of  photosynthesis; 
mol  carbon  m"3  s'1 

ratio  of  TCPY  to  maximum  temperature-adjusted  carbon  growth  yield 
cosine  of  the  average  downwelling  irradiance  angle 
cosine  of  the  average  upwelling  irradiance  angle 
Michaelis-Menten  half-saturation  constant  for  photon  saturation 
Attenuation  function  for  downwelling  irradiance,  direct  component;  m"1 
Attenuation  function  for  downwelling  irradiance,  diffuse  component;  m"1 
Attenuation  function  for  upwelling  irradiance;  m"1 
Irradiance  reflectance 
depth  level;  -meters 

absorption  fraction  per  depth  increment  (as  defined  in  the  text);  m"1 


Table  A3.  Optical  Subscripts 

sed  -  mineral  suspended  sediments 
det  -  particulate  organic  detritus 
cdm  -  colored  detrital  matter 
cdmb  -  background  cdm  signal 
cdom  -  colored  dissolved  organic  matter 
ph  -  phytoplankton 
p  -  particulates 

dir  -  direct  component  of  the  spectral  irradiance 
diff  —  the  diffuse  component  of  the  spectral  irradiance 
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